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We present a range-separated linear-response time-dependent density-functional theory (TDDFT) 
which combines a density-functional approximation for the short-range response kernel and a 
frequency-dependent second-order Bethe-Salpeter approximation for the long-range response kernel. 
This approach goes beyond the adiabatic approximation usually used in linear-response TDDFT and 
aims at improving the accuracy of calculations of electronic excitation energies of molecular systems. 
A detailed derivation of the frequency-dependent second-order Bethe-Salpeter correlation kernel is 
given using many-body Green-function theory. Preliminary tests of this range-separated TDDFT 
method are presented for the calculation of excitation energies of the He and Be atoms and small 
molecules (H 2 , N 2 , CO 2 , H 2 CO, and C 2 H 4 ). The results suggest that the addition of the long-range 
second-order Bethe-Salpeter correlation kernel overall slightly improves the excitation energies. 


I. INTRODUCTION 

Linear-response time-dependent density-functional 
theory (TDDFT) [1, 2] is nowadays one of the most 
popular approaches for calculating excitation energies 
and other response properties of electronic systems. 
Within the usual adiabatic semilocal density-functional 
approximations (DFAs), linear-response TDDFT usually 
provides reasonably accurate low-lying valence elec¬ 
tronic excitation energies of molecular systems at a low 
computational cost. However, these usual adiabatic 
semilocal DFAs have serious failures. In particular, 
they give largely underestimated Rydberg [3] and 
charge-transfer [4] excitation energies and they do not 
account for double (or multiple) excitations [5]. 

The problem with Rydberg and charge-transfer ex¬ 
citation energies is alleviated with the use of hybrid 
approximations in linear-response TDDFT [6], which 
combine a Hartree-Fock (HF) exchange response kernel 
with a DFA exchange-correlation response kernel. This 
problem is essentially solved with range-separated hy¬ 
brid approximations [7-10], introducing a long-range HF 
exchange kernel. Research in linear-response TDDFT 
now aims at an increasingly higher accuracy and re¬ 
liability, and in particular the inclusion of the effects 
of the double excitations. Examples of recent develop¬ 
ments are: the dressed TDDFT approach (combining 
TDDFT and the polarization-propagator approach) [11- 
13], double-hybrid TDDFT methods (combining TDDFT 
and configuration-interaction singles with doubles cor¬ 
rection [CIS(D)]) [14], and range-separated TDDFT ap¬ 
proaches in which the long-range response is treated with 
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density-matrix functional theory (DMFT) [15], multicon¬ 
figuration self-consistent-field (MCSCF) theory [16], or 
the second-order polarization-propagator approximation 
(SOPPA) [17]. 

In condensed-matter physics, the Bethe-Salpeter equa¬ 
tion (BSE) applied within the GW approximation (see, 
e.g., Refs. 18-20) is often considered as the most success¬ 
ful approach to overcome the limitations of TDDFT. Al¬ 
though it has been often used to describe excitons (bound 
electron-hole pair) in periodic systems, it is also increas¬ 
ingly applied to calculations of electronic excitation en¬ 
ergies in finite molecular systems [21-49]. In particular, 
the BSE approach was shown to give accurate charge- 
transfer excitation energies in molecules [30, 32, 34- 
36, 38, 39], and when used with a frequency-dependent 
kernel it is in principle capable of describing double exci¬ 
tations [33, 50, 51]. The drawbacks of the standard BSE 
approach with the usual approximations are the need to 
first perform a computationally demanding GW quasi¬ 
particle calculation, and an observed loss of accuracy for 
small molecules [44] which is probably due to self screen¬ 
ing. 

In this work, we explore the combination of TDDFT 
and BSE approaches based on a range separation of 
the electron-electron interaction. More specifically, we 
propose a range-separated TDDFT approach in which 
the long-range response is treated with a frequency- 
dependent second-order Bethe-Salpeter-equation (BSE2) 
correlation kernel. The BSE2 approximation was re¬ 
cently introduced by Zhang et al. [52] within the Tamm- 
Dancoff approximation (TDA) [53]. Compared to the 
standard BSE approach with the GW approximation, 
the BSE2 approximation keeps only second-order terms 
with respect to the electron-electron interaction, includ¬ 
ing second-order exchange terms which makes it free from 
self screening. It is an appropriate approximation for fi¬ 
nite molecular systems with relatively large gaps. Build¬ 
ing on the work of Sangalli et al. [51], we provide an alter¬ 
native and more general derivation of the BSE2 approxi¬ 
mation and we apply it to the range-separated case. We 
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present preliminary tests of this range-separated TDDFT 
method for the calculation of excitation energies of the 
He and Be atoms and some small molecules (H 2 , N 2 , 
C0 2 , H 2 CO, and C 2 H 4 ). 

In this range-separated TDDFT approach, an adia¬ 
batic semilocal DFA is used only for the short-range part 
of the exchange-correlation kernel, while a frequency de¬ 
pendence is introduced in the long-range part of the cor¬ 
relation kernel. This is motivated by the fact that the 
exact exchange kernel becomes spatially local and fre¬ 
quency independent in the limit of a very short-range 
interaction [10, 54], so that the adiabatic local-density 
approximation (LDA) becomes exact in this limit. Simi¬ 
larly, the short-range part of the exact correlation kernel 
is expected to be more spatially local and less frequency 
dependent than its long-range counterpart, so that an 
adiabatic semilocal DFA is expected to be accurate when 
restricted to the short-range part of the correlation ker¬ 
nel, as it happens for the ground-state correlation density 
functional [55]. 

Similarly to the ground-state case where second-order 
perturbation theory is appropriate for describing the 
long-range part of the correlation energy of systems with 
large enough gaps [56], the BSE2 approximation is ex¬ 
pected to be appropriate for describing the long-range 
part of the response of such systems. Moreover, in com¬ 
parison to the original full-range BSE2 scheme, the re¬ 
striction of the BSE2 approximation to the long-range 
part leads to potential practical and computational ad¬ 
vantages: (1) eliminating the need to do a first GW quasi¬ 
particle calculation since range-separated hybrid approx¬ 
imations provide orbital energies that are already close 
to quasiparticle energies [57, 58]; and (2) speeding up 
the computation of the BSE2 correlation kernel by us¬ 
ing multipolc expansions for the long-range two-electron 
integrals [59]. 

This paper is organized as follows. In Section II, 
we summarize the main equations of linear-response 
TDDFT with range separation. In Section III, we provide 
a full derivation of the frequency-dependent BSE2 corre¬ 
lation kernel without using the TDA, giving expressions 
in terms of space-spin coordinates and in a spin-orbital 
basis. Section IV explains how we practically perform 
the calculations and gives computational details for the 
systems tested. The results are given and discussed in 
Section V. Finally, Section VI contains our conclusions. 


II. RANGE-SEPARATED TIME-DEPENDENT 
DENSITY-FUNCTIONAL THEORY 

As a relatively straightforward extension of linear- 
response TDDFT [1], in range-separated TDDFT [10, 
16], the inverse of the frequency-dependent linear- 
response function is expressed as 

X _1 (xi,x 2 ;x , 1 ,X2;o;) = (x lr ) _1 ( x i, x 2 ; x' 1; x 2 ; w) 

~/hxc( X 1 ) x 2 ; x i , x 2 ; w), (1) 


where x = (r, a) stands for space-spin coordinates. 
In this expression, x' r ( x i, X 2 ; x) , x 2 ; w) is the linear- 
response function associated with the long-range (lr) in¬ 
teracting Hamiltonian 

H ll = f+V ne + Wll + V^ c , (2) 

where T is the kinetic-energy operator, Vne is 
the nuclei-electron interaction operator, W^ e is 
a long-range electron-electron interaction opera¬ 
tor, and Vg xc is a corresponding short-range (sr) 
Hartree-exchange-correlation (Hxc) potential op¬ 
erator. Additionally, /h xc ( Xi ’ x 2! x i) x 2 i cj) = 

/hxc(x1)X2; w)<5(xi,x , 1 )(5(x 2 ,X2) is the short-range 
Hxc kernel related to the functional derivative of the 
short-range Hxc potential with respect to the density 
(and 6 is the delta function). In practice, the long-range 
electron-electron interaction is defined with the error 
function as tyk(ri,r 2 ) = erf(/i|ri - r 2 |)/|ri - r 2 |, 
where the parameter /r controls the range of the 
interaction. Even though Eq. (1) is written with 
functions depending on four space-spin coordinates 
for generality, range-separated TDDFT only gives 
exactly the diagonal part of the linear-response function 
x(xi,x 2 ;w) = x( x i,x 2 ;xi,x 2 ;w), just as in usual 
TDDFT. 

In the time-dependent range-separated hybrid 
(TDRSH) scheme [10], the long-range linear-response 
function x' r ( w ) is calculated at the HF level. More 
precisely, the inverse of the long-range linear-response 
function is approximated as 

(x* r )~ 1 ( x i,x 2 ; x' 1 ,x' 2 ;oj) « (xo) _ 1 (xi,x 2 ;x , 1 ,x , 2 ;w) 

— /hx,hf ( x i 1 x 2 j x^ , x 2 ), (3) 

where Xo( w ) i s the non-interacting linear-response func¬ 
tion associated with the range-separated-hybrid (RSH) 
reference Hamiltonian [56] 

Ho = T + Vne + Vhx, HF + Vhxc, (4) 

with the long-range HF potential operator V^ r xH F> and 
/fj x ( x i) x 2 ! x i 5 x 2 ) is the corresponding long-range HF 
kernel. The latter is the sum of a long-range Hartree 
kernel, 

/h ( x i , x 2 ; x' x , x 2 ) = w(, r e (r 1 ,r 2 )5(x 1 ,x , 1 )5(x 2 ,x^), (5) 

and a long-range HF exchange kernel, 

fyL,HF ( X 1 1 x 21 X 11 x 2 ) = -^ r e( r l,r 2 )5(xi,X , 2)(5(x2,x' 1 ). (6) 

To go beyond the HF level, it was proposed to cal¬ 
culate x* r ( w ) a t the linear-response MCSCF level [16] 
or at the SOPPA level [17]. In the present work, we 
explore the recently proposed BSE2 approximation [52]. 
We thus propose to approximate the inverse of the long- 
range linear-response function as 

(x lr ) '( x i, x 2 ; x i, x 2 ;^) ~ (xo) -1 ( x i, x 2 ; x i, x 2 ;w) 

./hx.HF ( X 1; x 2 ! x l; x 2 ) 0c^BSE2 ( X 1 1 x 2 j x l) x 2 i ^0; 00 
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S< 2) (1,2) 



+ 



The Feynman diagrams of these terms are represented in 
Figure 1. 


B. Second-order Bethe-Salpeter correlation kernel 
in the time domain 


Figure 1: Feynman diagrams of the second-order correlation 
self-energy E( 2 \l,2). The time axis is vertical. The dots 
represent the outer variables 1 and 2. Horizontal dashed lines 
represent electron-electron interactions w ee - Arrowed lines 
represent one-particle Green functions G. The first diagram 
is the direct contribution of Eq. (9) and the second diagram 
is the exchange contribution of Eq. (10). 


with the long-range BSE2 frequency-dependent correla¬ 
tion kernel /c I BSE 2 ( a; ) f° r which we offer an alternative 
and more general derivation compared to Ref. 52. 


III. SECOND-ORDER BETHE-SALPETER 
CORRELATION KERNEL 

In this Section, we provide a derivation of the BSE2 
correlation kernel. For more details, see Ref. 60. We 
consider an arbitrary electron-electron interaction w ee in 
the derivation instead of the long-range one. 


A. Second-order correlation self-energy 


The starting point is the second-order correlation self- 
energy as a functional of the one-electron Green function 
G(l,2) where 1 = (xi,fi) and 2 = (x 2 ,f 2 ) stand for 
space-spin-time coordinates (see, e.g., Ref. 37) 


E c (2) (1,2)=* j d3d3'd4d4'd5d5' G(3,3') ^ 

tSee(3\ 4; 2, 4 ')xip( 4', 5; 4,5> e e(5', 1; 5,3), 

where w ee (l,2-,l',2') = izj ee (l, 2)5(1, l')5(2, 2') is an ar¬ 
bitrary electron-electron interaction, w ee (l, 2; 1', 2') = 
u; ee (l, 2; 1', 2') — ui ee (2,1; 1', 2') is the corresponding 
antisymmetrized interaction, and xip(l, 2; 1', 2') = 

—*G(1, 2')G(2,1') is the independent-particle (IP) four- 
point linear-response function. The presence of the an¬ 
tisymmetrized interaction w ee in Eq. (8) means that the 
second-order correlation self-energy can be decomposed 
as Ec 2 '* = Ec 2d ^ + Ec 2x) with a direct contribution 


E< 2d )(l,2) = *G( 1,2) y’d3d4«)ee(2,3) 
x Xip( 3, 4; 3,4)uj ee (4,1), 
and an exchange contribution 

E( 2x )(1,2) = -i J d3d4 G(l, 3)u> ee (2, 3) 
x Xip(3, 4; 2, 4)iu e e(4, l). 


(9) 


The second-order Bethe-Salpeter correlation kernel is 
defined as the functional derivative of the second-order 
correlation self-energy with respect to the Green function 


S( 2) (l,4;2,3) 


.<5Ec 2) (1,2) 
* JG(3,4) 


( 11 ) 


Taking the derivative of Eq. (8) generates three terms 


S< 2) (1,4; 2,3) = 

- J d5d5'd6d6'u)ee(4,5; 2, 5 ')xip(5', 6; 5,6')ui e e(6', 1; 6,3) 

- J d5d5'd6d6'u)ee(5,4; 2, 6 )xip( 5', 6; 6', 5)w ee (6', 1; 3,5') 

- f d5d5'd6d6'w; ee (5,6; 2, 3 )xip( 6', 5'; 6, 5)w; ee (4,1; 5', 6'), 

( 12 ) 


which, as done for the correlation self-energy, can be de¬ 
composed as Sc 2 ' 1 = Sc 2d) + si 2x -*, with a direct contribu¬ 
tion 

£( 2d >(l,4;2,3) = 

- <5(1,3)<5(2,4) J d5d6 w ee (2, 5)xip( 5, 6; 5,6)w ee (6,1) 

- w ee ( 2, 4)xip( 1, 4; 3, 2)w ee (3,1) 

- w ee ( 2, 3)xip( 1, 4; 3, 2) Wee (4,1), (13) 

and an exchange contribution 

£( 2x )(l,4;2,3) = 

5(1,3) J d5ta ee (2,4) X ip(4,5; 2,5)w ee (5,1) 

+ 5(2,4) j d5w ee (2, 5 )xip( 1, 5; 3, 5)m ee (3,1) 

+ u; e e(2,3)xip(l,4;2,3Ke(4,l). (14) 


The Feynman diagrams of these six terms are represented 
in Figure 2. Similar kernel diagrams are shown in Ref. 61. 

Introducing explicitly the time variables, using an 
instantaneous spin-independent electron-electron in¬ 
teraction in ee (l,2) = u> ee (iT, r 2 )5(fi, t 2 ) and time- 
translation invariance, we found that the second-order 
Bethe-Salpeter correlation kernel is composed of a 
particle-hole/hole-particle (ph/hp) part and a particle- 
particle/hole-hole (pp/hh) part, which non-trivially de¬ 
pend on only one time difference t\ — t 2 , 

5^ 2) (xiti,X4t4;X2t2,X3t 3 ) 

= 5(ti, t 3 )S(t 2 , f 4 )S( 2 ’ pll/hp) (xi, x 4 ; x 2 , x 3 ; ti - t 2 ) 

+ 5(ti, t 4 )5(f 2 , t 3 )5(. 2 ’ pp/hh) (xi, x 4 ; x 2 , x 3 ; ti - t 2 ), 

(15) 


( 10 ) 
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Figure 2: Feynman diagrams of the second-order Bethe-Salpeter correlation kernel Ec 2 ^(l, 4; 2, 3). The three upper diagrams 
are the direct contributions of Eq. (13) and the three lower diagrams are the exchange contributions of Eq. (14). The four 
diagrams on the left correspond to ph/hp terms and the two diagrams on the right correspond to pp/hh diagrams. Since the 
kernel is the functional derivative of the self-energy with respect to the Green function, these diagrams can be obtained from 
the ones of Figure 1 by removing one of the arrowed lines. 


with the ph/hp kernel 

s (2, ph/hp) ( Xl)X4;X2)X3;T ) = 

= - (5(xi,x 3 )<5(x2,x 4 )/dx 5 d We (r 2 ,r 5 ) 

x Xip( x 5,x 6 ;x 5 ,x 6 ; -r)uJee(r6,ri) 

- w ee (r 2 , r 4 )xip (xi, x 4 ; x 3 , x 2 ; T)w ee (r 3 , ri) 

f (16) 

+ 5(xi,x 3 ) / dx 5 u>ee(r 2 ,r 4 ) 

x xip(x 4 ,x 5 ;x2,x 5 ;-r)uJee(r5,ri) 

+ J(x 2 ,x 4 ) /dWr 2 ,r 5 ) 
x Xip( x 1 )x 5 ;x 3 ,x 5 ;r)wee(r 3 ,r 4 ), 


and the pp/hh kernel 
7 -( 2 pp/hh) (x 4 , x 4 ; x 2 , x 3 ; t) 

= - Wee(r 2 , r 3 )Xip /hh ( x i> x 4; x 3, x 2 ; T)u; ee (r 4 , r 4 ) (17) 
+ Wee (r 2 , r 3 )xl > p /llh (x 4 , x 4 ; x 2 , x 3 ; r)w ee (r 4 , ri), 


domain which is [62, 63] 


x(1,2;1' j 2')=Xip(1,2;1' ) 2') + J d3d4d5d6 

Xip(1, 4; T, 3)S H xc( 3, 6; 4, 5)x(5, 2; 6, 2'), 

(18) 

where x(l,2;l',2') is the interacting four-point linear- 
response function and Sjjxc is the Bethe-Salpeter Hxc 
kernel. Written explicitly with time variables, and setting 
t\ = tf and t' 2 = t 2 (where t + = t + 0 + refers to a time 
variable with an infinitesimal positive shift) to extract the 
(ph/hp) linear-response function, the equation becomes 

X(x 4 ti, x 2 t 2 ; , x' 2 t^) = Xip( x iii> x 2 i 2 ; x 'i^, x 2 ^) 

+ J dx 3 dt 3 dx 4 dt 4 dx 5 dt 5 dx 6 dt 6 Xip( x iii, x 4 t 4 ; x] tf, x 3 t 3 ) 

S H xc( x 3i3, x 6^6; x 4<4, X5t 5 )x(x 5 t 5 , X 2 t 2 ; X 6 t 6 , X^), 

(19) 


where xip ( x i 7 x 2 5 x 'i, x 2 ; t = ti - t 2 ) = 

Xip(xit 4 ,x 2 t 2 ;x/ti,x 2 t 2 ) is the IP (ph/hp) linear- 

response function and Xip hh ( x i,x 2 ; x' 1; x 2 ; r = t 4 — 1 2 ) = 
Xip(x 4 ti, x 2 t 4 ; x.[t 2 , x 2 t 2 ) is the IP pp/hh linear-response 
function. As the names suggest, the former describes the 
independent propagation of one particle and one hole, 
and the latter describes the independent propagation 
of either two particles or two holes, depending on the 
sign of t\ — t’ 2 - Because of the different delta functions 
on the time variables in Eq. (15), the ph/hp and pp/hh 
contributions need to be treated separately. 

C. Effective second-order Bethe-Salpeter 
correlation kernel in the frequency domain 

The Bethe-Salpeter kernel that we derived must be 
used in the general Bethe-Salpeter equation in the time 


where S Hx c( x 3 f 3 , x-sU] x 4 ^ 4 , x 5 ^) = 

(o\ 

/hx.hf ( x 3 1 x 6; x 4 1 x 5 ) + ; (x 3 t 3 ,x 6 t 6 ;x 4 t 4 ,x 5 t 5 ) is 

taken as the sum of time-independent HF kernel /hx.hf 

and the second-order correlation kernel 3c. Because 
of the time dependence in Sc , in Eq. (19) the time 
variables t 3 and f 4 cannot be equated, and neither can 
be the time variables £5 and t§. Consequently, Eq. (19) 
is not a closed equation for the (ph/hp) linear-response 
function x( x iii, x 2 i 2 ; x 2 i 2 )■ 

To close the equation, Zhang et al. [52] followed Stri- 
nati [18] and used an explicit time-dependent form in 
the TDA for the ph/hp amplitudes [18, 64] with which 
X in Eq. (19) can be expressed. Here, instead, following 
Sangalli et al. [51] (see also Ref. 50), we work in Fourier 
space and define an effective kernel depending on only 
one frequency without using the TDA. Introducing the 
decomposition of 3/ ; in ph/hp and pp/hh terms given 
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in Eq. (15) and Fourier transforming Eq. (19) gives 
xM = XipM + Xip M /hx xM 

+ J ^r^-xip ( W ', W )S( 2 ’P h / h P)( W '- W ")x(^^) 

+ J ^ ^Xip( u,»~( 2 ' pp / hh ) (u/ + /) X («/, U ), 

( 20 ) 

where the space-spin variables have been dropped for 
conciseness (all the quantities depend on four space-spin 
variables), and the integrations over uj' and u>" are from 
—oo to +oo. In this expression, x(w , ,w) is the double 
Fourier transform 

X(w',w) = J dndr e i ^e i ^x(n,T 2 = 0 ~,r), (21) 


where x( r i> T 2 , t) corresponds to the response function 
x(x-] t\ , x 2 f 2 ; x( t\ , x 2 f 2 ) expressed with the time variables 
ri = T 2 = t 2 -t' 2 , and r = (ti+t , 1 )/2-(t 2 +t 2 )/2, 

and similarly for xip(u/, ui). As a special case, x(w) is 
just the Fourier transform of the (ph/hp) linear-response 
function x( T i = 0, T 2 = 0, r), and similarly for xip( w )- 
Obviously, ~;(. 2 ’ ph / hp ) ( w ) an d w( 2 > pp / hh ) ( w ) are the Fourier 
transforms of Sc 2,ph ^ hp '(r) and ~;( 2 ’ pp / hh ) ( T ) given in 
Eqs. (16) and (17), respectively. Eq. (20) can be rewrit¬ 
ten as an effective Bethe-Salpeter equation involving only 
one frequency [51] 


xM =xipM + xipM/hxxM + xipM4 2) MxM, 

( 22 ) 

or, equivalently, 

X~V) = XipV) - /hx - 2( 2) (w), (23) 

with an effective correlation kernel defined as 

do/ d u" 


s i 2) M = Xip 1 ( w ) 
n( 

1 —'c 

+ Xip 1 M 


, 2n 2, x,p ^^ 
S(>,ph,bp) ( „- -„» )x( „",„ )x -'(„ ) 

dw'dw" . , , 

-Xip(w ,w) 


27T 27T 

5( 2 ’ pp/hh) (w , + w")x(w ,/ ,w)x" 1 (w). 


(24) 


To keep only second-order terms in Eq. (24) we must re¬ 
place both the IP linear-response function xip and inter¬ 
acting linear-response function x by the non-interacting 
linear-response function xo ; and we finally arrive at the 
BSE2 correlation kernel 


/c,bse 2 M = XoV) J ^ ^ Xq(u',w) 

s (2,ph/hp) (w ,_^ )xo(u; » ;W)x -l (w) 

_i , . /■ du/dw" , 

+X ° ^r^ xo(w ’ w) 

2(2.P p /hh) (w / + w //) Xo ( U ;", w )x- 1 (a;), 


( 25 ) 


where ~f 2,ph / hp ) and ;=;( 2 > pp / hh ) are obtained from 
Eqs. (16) and (17) with the replacement of xip by xo 
as well. 

We note that, in Eq. (23), xip 1 (w) could also be ex¬ 
panded to second order, leading to self-energy (or quasi¬ 
particle) contributions to the effective kernel [51]. How¬ 
ever, in this work, we do not consider such self-energy 
contributions to the kernel. 


D. Expressions in a spin-orbital basis 


We now give expressions in the orthonormal canonical 
spin-orbital basis {<p P } of the reference non-interacting 
Hamiltonian. Any function i r, (xi,x 2 ;x' 1 ,x 2 ) depending 
on four space-spin coordinates can be expanded in the 
basis of products of two spin orbitals, and its matrix el¬ 
ements are defined as 


1 7)1 



i dx' x dx 2 dx 2 (x'i ) ip* (xi ) 


F (xi, x 2 ; xi, x!, ) ip* (x 2 ) <p s (x !,), 


(26) 


where p, q , r, s refer to any (occupied or virtual) spin or¬ 
bital. In the following, the indices i,j,k,l will refer to 
occupied spin orbitals and the indices a , b, c, d to virtual 
spin orbitals. 

Using the expression of the Fourier transform of the 
non-interacting (ph/ph) linear-response function, 


Xo(xi,x 2 ;xi,x' 2 ;w) = 

y- Vl (x'i ) ‘fc (xi) ip* (x' 2 ) ip k (x 2 ) 
kc <*> - (e c - e fc ) + *0+ (27) 

\ - ( / ? fc(x 2 )y ; c(x 2 )^* (x / 1 )l/?fc(xi) 

“ uj + (e c - e k ) - i0 + 

and of the non-interacting pp/hh linear-response func¬ 
tion, 


^.PP/hh 

Xo 


(xi,x 2 ;xi,x 2 ;w) = 


E 


y j fc( x 'i)^( x i)y ? r( x 2)yfc(x 2 ) 

w — (£fc + e l) ~ *0+ 


-E 


y j c( x 'i)^( x i)^( x 2)y ; c(x2) 

w + (e c + £ d) + *0+ 


(28) 


where e p are the spin-orbital energies, we find the matrix 
elements of the Fourier transform of the ph/hp second- 
order correlation kernel, 

=l2.ph/hpW.o _ V- (rc\\pk)(kq\\cs) 

A™-™ <'' , -A«-( e .-e 1 )+# 

X- efc||pc){cq||fc5) 

“ u + (e c - Ek) - i0+ ’ 
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and of the pp/hh second-order correlation kernel, 


the linear-response matrix usually denoted by B), we find 


-( 2 ,pp/hh)/ 1 y- (qr\\kl){lk\\sp) 

-c, pq , rs 2 Z^ w _ (£fc+£;) _ z0+ 

1 y (qr\\cd)(dc\\sp) 

2 w _ (e c -p e d ) + *0+ 


(30) 


where (pg||rs) = (pq\rs) — (pq\sr) are the antisym¬ 
metrized two-electron integrals associated with the in¬ 
teraction w ee . 

The matrices of xo(w) and yo(w / , ui) are both diagonal 
with elements in the occupied-virtual/occupied-virtual 
spin-orbital product block given by 


Xo,; 


t (w) = 


- (e Q - Si) + i0+ ’ 


(31) 


and (see Appendix A) 


X0,ia,ia(bJ ;^) — i 6 XO 

1 1 


i T ui /2 — s a T i0^~ ui f — u) j 2 — Ei — 


,(32) 


and, for the virtual-occupied/virtual-occupied block, 

Xo ,ai,ai(w',U!) = X0,m,ia(w', -U>) and X'O,ai,oi(w) = 

Xopa,ia( —1 w). The matrix elements of the BSE2 corre¬ 
lation kernel are then found straightforwardly by do¬ 
ing the matrix multiplications and contour-integrating 
over the frequencies in the upper-half complex plane 
in Eq. (25). For the matrix elements in the occupied- 
virtual/occupied-virtual (ov/ov) block (contributing to 
the linear-response matrix usually denoted by A), we 
find 


/c,BSE2,ia, t&M 


_ (jc\\ik)(ka\\cb) 

“ W - (£ b + 6 C - £i - e k ) 

_ y- (jk\\ic)(ca\\kb) 

“ w - (e a + e c - £j - £fe) 

+ 1 (aj\\kl)(lk\\bi) 

2 ^ ui - (e a + E b - £k - £i) 

| 1 {aj\\cd)(dc\\bi) 

2 ■— w — (e c + £d — £i — £?) ’ 


(33) 


fc,BSE2,ia,bj ^ ) 

kc 

-E 


fee 


(bc||ifc)(fca||cj) 

— (Eb + E c — Ei — Ek) 

(bk\\ic)(ca\\kj) 

— (e a + e c — — £fe) 


| 1 (a,6||fcZ)<ZA||jz) 

2 y-' -(e a +£& -£fc -£j) 


| 1 (ab\\cd)(dc\\ji) 

2^ -{Ec + Ed-Ei-EjY 


(34) 


which turn out to be independent of the frequency. 
To the best of our knowledge, the matrix elements in 
Eq. (34) had never been given in the literature before. 
It is easy to check that the ov/ov block is Hermitian, 
/c,BSE 2 ,iajb(w) = / c ,BSE2, jb,ia M*, and that the ov/vo 
block is symmetric, f c ,BSE2,ia,bj = fc,BSE2,jb,ai- 

The matrix elements of the BSE2 correlation kernel 
display sums over either one occupied and one virtual or¬ 
bital (for the ph/hp terms) or over two occupied or two 
virtual orbitals (for the pp/hh terms). In a straightfor¬ 
ward implementation, the computational cost of the lat¬ 
ter scales as A/ Ay) where N a is the number of occupied 
orbitals and N v the number of virtual ones. However, 
in the case of the long-range interaction, the computa¬ 
tional cost of the BSE2 correlation kernel could be made 
low, e.g. by approximating the long-range two-electron 
integrals by multipole expansions [59]. 


IV. PRACTICAL RESOLUTION AND 
COMPUTATIONAL DETAILS 

A. Perturbative resolution 

In the range-separated scheme that we propose, we ap¬ 
proximate the inverse of the linear-response function as 
[combining Eqs. (1) and (7)] 

X : (w) « Xo X ( w ) — /hx,hf - /hxc - /c)bse 2( w ), (35) 


Note that the denominators of Eq. (33) contain the sum 
of two virtual spin-orbital energies minus the sum of 
two occupied spin-orbital energies, i.e. a non-interacting 
double-excitation energy. Thus, the denominators are 
small (and therefore the kernel can be large) whenever 
ui is close to a non-interacting double-excitation energy. 
The matrix elements in Eq. (33) are identical (at least for 
real-valued spin orbitals) to the kernel matrix elements 
recently derived by Zhang et al. [52] in the TDA [65]. 
The matrix elements in Eq. (33) also show some simil¬ 
itude with the SOPPA kernel [12, 66-68] and the sec¬ 
ond RPA kernel [51]. Similarly, for the matrix ele¬ 
ments of the BSE2 correlation kernel in the occupied- 
virtual/virtual-occupied (ov/vo) block (contributing to 


where Xo( w ) the R.SH non-interacting linear-response 
function and // r bse 2 ( w ) is the BSE2 correlation kernel 
for the long-range electron-electron interaction. We note 
that, according to Eq. (23), instead of Xo 1 (w), we should 
use in Eq. (35) the inverse of the long-range IP linear- 
response function (x}p) _1 ( w ) constructed with the long- 
range interacting Green function. This could be ac¬ 
counted for by either adding quasiparticle corrections to 
the orbital energies, as done in Ref. 52, or adding self¬ 
energy contributions to the long-range correlation ker¬ 
nel [51]. These contributions can generally be important 
when using HF orbitals or DFT orbitals with semilocal 
DFAs. However, in the case of range separation, the 
orbital energies obtained with long-range HF exchange 
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are already good approximations to quasiparticle ener¬ 
gies [57, 58]. It is thus reasonable to use the approxima¬ 
tion (xjp) -1 (w) « Xq 1 (w). We come back to the possi¬ 
bility of adding quasiparticle corrections in Section and 
discuss their effects on He, Be, and H 2 in Section V A. 

When projected in the basis of the RSH spin orbitals, 
Eq. (35) leads to the self-consistent pseudo-Hermitian 
eigenvalue equation 

( AK) B \(X n \_ fl 0 \fX n \ 

^ B* A(-w„)* ) ^ Y„ ) ~ W " V 0 -1 )\Y n )' 

(36) 

where uj n are the excitation (or diexcitation) energies, 
(X„,Y„) are the associated linear-response eigenvectors, 
and the matrix elements of A and B are given by 

A ia j b (uj) = (e a - £i)5ij5 ab + ( aj\w ee \ib ) - (aj\w^ e \bi) 

+fxc,ia,jb + fc,BSE2,ia,jb( UJ ): ( 37 ) 

and 

B ia ,jb = ( ab\w ee \ij } - (ab\w l g e \ji) 

+ fxc,ia,bj + /c!bSE2, ia,bji (38) 

where e p are the RSH spin-orbital energies, (pq\wee\rs) 
and (pq\w l * e \rs) are two-electron integrals in the RSH 
spin-orbital basis associated with the Coulomb interac¬ 
tion w ee and the long-range interaction ui* r e , respectively, 
and /xc lP9i r S are the matrix elements of the short-range 
exchange-correlation kernel. The matrix elements of the 
long-range BSE2 correlation kernel /c][bse 2 ,pg,rs are gi ven 
in Eqs. (33) and (34) using in these expressions long- 
range two-electron integrals (pq\\rs) —> (pq|?z;*g|rs) — 
(pg|w;g r e |sr) and RSH spin-orbital energies e p . 

The resolution of the self-consistent eigenvalue equa¬ 
tion (36) is more complicated than in the standard case 
of a frequency-independent matrix A. Following Zhang 
et al. [52], for a first exploration of the method, we work 
within the TDA (i.e., we set B = 0) and use a non-self- 
consistent perturbative resolution. We thus decompose 
the matrix A in Eq. (37) as the sum of the frequency- 
independent RSH contribution [10] and the long-range 
frequency-dependent BSE2 correlation kernel contribu¬ 
tion 

A (uj) = ArsH + fc(BSE2( w )- (39) 

The TDRSH linear-response equation is first solved in 
the TDA, 

ArshXo,^ = wo, n X 0jn , (40) 

where cj o, n and Xo >n are the corresponding excitation en¬ 
ergies and linear-response eigenvectors, respectively. The 
effect of the long-range BSE2 correlation kernel is then 
added perturbatively to obtain the excitation energies 

= Cdo,n + z n xj n fc]BSE 2 ( a; 0 ,n) Xo jK , (41) 


where Z n is the normalization factor 


II 

1 ^c(bSE 2( w ) 

1 X °’" dco 

XoJ 


V 

^~^0,n ) 


As pointed out by Zhang et al. [52], the effect of the nor¬ 
malization factor Z n turns out to be very small (Z n is 
always very close to 1, especially in the range-separated 
case), but we keep it in our calculations. We note that 
the expression of the correction Xj n fg r BSE2 (wo,n) Xo,„ 
in Eq. (41) is very similar (but not identical) to the 
so-called “direct” contribution of the CIS(D) correc¬ 
tion [69, 70]. As for CIS(D), it is easy to check that 
Xo n fc r BSE 2 ( w o,n) Xo,n contains only connected terms 
and thus provides a size-consistent correction to the exci¬ 
tation energies. Using this non-self-consistent perturba¬ 
tive resolution has the consequence that the total num¬ 
ber of calculated excitation energies is equal to the num¬ 
ber of single excitations, so we cannot obtain excita¬ 
tions with primarily double-excitation character. How¬ 
ever, the BSE2 correlation kernel brings the effects of 
non-interacting double excitations on excited states with 
dominant single-excitation character. 

The method defined by Eqs. (40) and (41) will be re¬ 
ferred to as TDRSH+BSE2. When the range-separation 
parameter p is set to zero, all long-range contributions 
vanish, and it reduces to the standard time-dependent 
Kohn-Sham (TDKS) method in the TDA. When p goes 
to Too, all short-range contributions vanish, and it re¬ 
duces to time-dependent Hartree-Fock (TDHF) within 
the TDA [i.e., configuration-interaction singles (CIS)] 
with a BSE2 correction, which will be referred to as 
TDHFTBSE2. As regards the density-functional ap¬ 
proximation, in this work, we use the short-range LDA 
exchange-correlation functional of Ref. 71 in the ground- 
state RSH calculations (i.e., for determining the RSH or¬ 
bitals and orbital energies) and the corresponding short- 
range LDA exchange-correlation kernel [10] in the linear- 
response TDRSH calculations. 

B. Long-range excitation energies 

For He, Be, and H 2 , we also perform calculations of 
long-range excitation energies as a function of p (i.e., 
along the range-separated adiabatic connection, similarly 
to Refs. 72-74) obtained by removing the contribution 
from the short-range Hxc kernel / Bxc in the matrix ele¬ 
ments A ia j b of Eq. (37), i.e. 

= (e a - £i)6ij6 ab T (aj\w l J e \ib) - (aj\w l * e \bi) 

+fc,BSE2,ia,jb( UJ ), ( 43 ) 

within the perturbative resolution of Eqs. (40) and (41) 
in the TDA. The orbitals and orbital energies used in 
Eq. (43) are still the RSH ones (i.e., with the short-range 
LDA exchange-correlation functional), as for the other 




calculations. The obtained long-range excitation ener¬ 
gies are approximations to the excitation energies of the 
long-range interacting Hamiltonian of Eq. (2), which re¬ 
duces to the LDA orbital energy differences at /Li = 0 
and to the TDHF+BSE2 excitation energies for /x —> oo. 
These long-range excitation energies allows us to test the 
effect of the BSE2 correlation kernel independently of 
the approximation used for the short-range exchange- 
correlation kernel, since we have accurate reference values 
for these quantities from Ref. 72. 

For these systems, we also test the addition of the per¬ 
turbative quasiparticle correction using the long-range 
second-order correlation self-energy, similarly to Ref. 52, 
i.e. replacing the RSH orbital energies e p in Eq. (43), 
including in the long-range BSE2 correlation kernel 
/c"BSE2,ia,jb( w )> by the quasiparticle energies 

e p = e p + z p E CiPp ( £p ), (44) 

with the renormalization factor z p = [1 — 

(dEc^pM/dw)^^]- 1 . In Eq. (44), E(f pp (e p ) is 
the diagonal matrix element of the frequency-dependent 
long-range second-order correlation self-energy E(, r (w) 
over the RSH spin orbital <^> p (x) evaluated at u = e p , 
whose expression is 


yUr 
^c,pp 


(w) 


1 \(ab\w^ e \pi) - (ab\w^ e \ip)\ 2 

2 W + Ei - £ a - £ b 

iab 

+ l y-' Ifa'KeM - fa>ie|ap)| 2 

2 LU -(- £ a £i £j 

tja J 


where i. j and a, b refer to occupied and virtual RSH spin 
orbitals, respectively. This quasiparticle correction will 
be denoted by GW2 since it is a second-order GIF-type 
correction. The resulting method will thus be referred to 
as GW2+TDHF+BSE2. 


C. Computational details 

We calculate vertical excitation energies of four small 
molecules, N 2 , CO, H 2 CO, and C 2 H 4 , at their experimen¬ 
tal geometries [75-78], using the Sadlej+ basis sets [79]. 
Our reference values are obtained by equation-of-motion 
coupled-cluster singles doubles (EOM-CCSD) calcula¬ 
tions performed with GAUSSIAN 09 [80]. For each 
molecule, we report the first 14 excited states found with 
the EOM-CCSD method. For each molecule, we per¬ 
form a self-consistent ground-state RSH calculation using 
the short-range LDA exchange-correlation functional of 
Ref. 71, followed by a spin-adapted closed-shell TDRSH 
linear-response calculation in the TDA using the short- 
range LDA exchange-correlation kernel [10], as imple¬ 
mented in a development version of MOLPRO [81]. The 
TDRSH+BSE2 excitation energies are then calculated 
by a spin-adapted closed-shell version of Eq. (41) imple¬ 
mented in a homemade software interfaced with MOL¬ 
PRO (see Ref. 60 for details). The range-separation pa¬ 
rameter /x is set to 0.35 bohr -1 which yields a minimal 


mean absolute deviation (MAD) over the four molecules 
of the TDRSH+BSE2 excitation energies with respect 
to the EOM-CCSD references. We note that it has 
been proposed to adjust the value of /x for each system 
by imposing a self-consistent Koopmans’ theorem condi¬ 
tion [82, 83] or, equivalently, minimizing the deviation 
from the piecewise linearity behavior of the total energy 
as a function of the electron number [84, 85]. This ap¬ 
proach is appealing but it has the disadvantage of be¬ 
ing non size consistent [86], so we prefer to use a fixed 
value of /j, independent of the system. For comparison, 
we also perform standard, linear-response TDKS calcu¬ 
lations with the LDA functional [87], as well as TDHF 
and TDHF+BSE2 calculations, all in the TDA. In the 
TDA, Xq t n t i a can be considered as the coefficient of the 
(spin-orbital) single excitation i —> a in the wave func¬ 
tion of the excited state n. Each excited state was thus 
assigned by looking at its symmetry and at the leading 
orbital contributions to the excitation. 

The calculations of the long-range excitation ener¬ 
gies for He, Be, and H 2 are done similarly except that 
the short-range LDA exchange-correlation kernel is re¬ 
moved in the TDRSH linear-response calculation. The 
GW2 quasiparticle correction is calculated using a spin- 
adapted closed-shell version of Eq. (44). We use an un¬ 
contracted t-aug-cc-pV5Z basis set for He, an uncon¬ 
tracted d-aug-cc-pVDZ basis set for Be, and an uncon¬ 
tracted d-aug-cc-pVTZ basis set for H 2 , for which we 
have reference long-range excitation energies obtained at 
the full configuration-interaction (FCI) level using an ac¬ 
curate Lieb-optimized short-range potential [72]. 


V. RESULTS AND DISCUSSION 

A. Long-range excitation energies of the He and 
Be atoms and of the H 2 molecule 

The long-range excitation energies to the first triplet 
and singlet excited states of the He atom are plotted as 
a function of the range-separation parameter fj, in Fig¬ 
ure 3. The triplet and singlet excitation energies are iden¬ 
tical at /x = 0, where they reduce to the non-interacting 
Kohn-Sham excitation energies. When increasing /x, i.e. 
when adding the long-range interaction, this degeneracy 
is lifted and the excitation energies tend to the physical 
excitation energies in the limit /x — > oo. At /lx = 0, for 
all the approximate methods tested here, the long-range 
excitation energies reduce to LDA orbital energy differ¬ 
ences, which, as well known for Rydberg states, strongly 
underestimate the exact Kohn-Sham orbital energy dif¬ 
ferences (by about 5 eV in the present case). This under¬ 
estimation of the long-range excitation energies is pro¬ 
gressively eliminated by increasing the value of /x until 
/x « 1 bohr -1 . For /x > 1.5 bohr -1 , with all the ap¬ 
proximate methods, the long-range excitation energies 
vary much less and are a bit too high compared to the 
reference FCI long-range excitation energies. The BSE2 
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Figure 3: Long-range excitation energies to the first triplet 
and singlet excited states of the He atom as a function 
of the range-separation parameter /x, obtained by long- 
range TDHF, long-range TDHF+BSE2, and long-range 
GW2+TDHF+BSE2 calculations in the TDA using RSH 
(with the short-range LDA functional) orbitals and an un¬ 
contracted t-aug-cc-pV5Z basis set. The reference FCI long- 
range excitation energies are from Ref. 72. 



Range-separation parameter p (bohr' 1 ) 

Figure 4: Long-range excitation energies to the first triplet 
and singlet excited states of the Be atom as a function 
of the range-separation parameter /x, obtained by long- 
range TDHF, long-range TDHF+BSE2, and long-range 
GW2+TDHF+BSE2 calculations in the TDA using RSH 
(with the short-range LDA functional) orbitals and an un¬ 
contracted d-aug-cc-pVDZ basis set. The reference FCI long- 
range excitation energies are from Ref. 72. 


correlation kernel has almost no effect for the singlet ex¬ 
cited state, while it increases the excitation energy for 
the triplet excited state which leads to a larger error at 
large /x. The GW2 quasiparticle correction systemati¬ 
cally decreases the excitation energies, leading to smaller 
errors at large /x for both singlet and triplet excitation 
energies. The GW2 correction on the excitation energies 
is relatively large (0.5 eV) for large /x, but decreases when 
/x is decreased, being less than 0.2 eV for /x < 1 bohr -1 
and about 0.01 eV for /x = 0.35 bohr -1 (the value of /x 
used for the other systems in Section VB). 

Figure 4 shows the long-range excitation energies to 


Figure 5: Long-range excitation energies to the first triplet 
and singlet excited states of the H 2 molecule at the equi¬ 
librium internuclear distance as a function of the range- 
separation parameter /x, obtained by long-range TDHF, long- 
range TDHF+BSE2, and long-range GW2+TDHF+BSE2 
calculations in the TDA using RSH (with the short-range 
LDA functional) orbitals and an uncontracted d-aug-cc-pVTZ 
basis set. The reference FCI long-range excitation energies are 
from Ref. 72. 


the first triplet and singlet excited states of the Be atom. 
For these low-lying valence states, the TDHF long-range 
excitation energies are relatively accurate close to /x = 0, 
but they deteriorate somewhat as /x is increased. For 
/x > 1 bohr -1 , TDHF underestimates the triplet long- 
range excitation energy by about 1 eV and the singlet 
long-range excitation energy by about 0.25 eV, in com¬ 
parison to the reference FCI long-range excitation en¬ 
ergies. Adding the BSE2 correlation kernel correctly in¬ 
creases the TDHF long-range excitation energies, leading 
to a relatively accurate singlet excitation energy for all /x 
and reducing the error in the TDHF triplet excitation en¬ 
ergy by a factor of 3 for large /x. The GW2 quasiparticle 
correction does not change much the excitation energies, 
at most about 0.1 eV. 

Finally, the long-range excitation energies to the first 
triplet and singlet excited states of the H 2 molecule at 
the equilibrium internuclear distance are reported in Fig¬ 
ure 4. For these molecular valence states, the LDA or¬ 
bital energy differences at /x = 0 are too low by more 
than 1 eV. Again, this underestimation is corrected by 
increasing the value of /x. For /x > 1 bohr -1 , TDHF al¬ 
ways underestimates the triplet long-range excitation en¬ 
ergy, while it is more accurate for the singlet long-range 
excitation energy. The addition of the BSE2 correlation 
kernel changes little the excitation energy for the singlet 
state, and significantly reduces the error on the excitation 
energy for the triplet excited state. The GW2 quasipar¬ 
ticle correction tends to improve a bit the accuracy of 
the excitation energies for intermediate values of /x, but 
remains very small for small and large values of /x (in 
particular, it is about 0.05 eV for /i = 0.35 bohr -1 ). 

Overall, these results are encouraging and support the 
relevance of the TDHF+BSE2 approximation for the 
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Table I: Excitation energies of N 2 calculated by linear-response TDKS (with the LDA functional), TDRSH and TDRSH+BSE2 
(with the short-range LDA functional and fi = 0.35 bohr^ 1 ), TDHF and TDHF+BSE2, all within the TDA. The EOM-CCSD 
excitation energies are taken as reference. The Sadlej+ basis set is used. 


State 

Transition 

TDKS 

TDRSH 

TDRSH+BSE2 

TDHF TDHF+BSE2 

EOM-CCSD 

Valence excitation energies (eV) 


l7T u 

->• 17r g 

8.08 

7.74 

7.93 

6.23 

8.88 

7.72 

3 n g 

3<Tg 

— >• l7Tg 

7.58 

7.85 

8.05 

7.99 

10.97 

8.16 

3 a u 

l7T u 

-A l7r (S 

8.88 

8.54 

8.74 

7.32 

9.96 

9.07 

X 

3f7 g 

-A l7r g 

9.17 

9.50 

9.68 

10.02 

12.43 

9.55 

3 W 

l7T u 

-A l7r g 

9.65 

9.34 

9.53 

8.50 

10.77 

10.00 

Wf 

l7T u 

-A l7r g 

9.65 

9.34 

9.53 

8.50 

10.84 

10.24 

3 a u 

l7T u 

-A l7r g 

10.25 

9.98 

10.18 

9.06 

11.30 

10.66 

3 n u 

2o-u 

-A l7r K 

10.42 

10.77 

10.97 

11.74 

14.82 

11.36 

Rydberg excitation energies (eV) 

3 

3f7 g 

—S- 4cr g 

10.28 

11.47 

11.56 

13.12 

13.94 

11.74 

3f7 g 

—> 4cr g 

10.40 

11.94 

11.98 

14.01 

14.22 

12.15 

3<r g 

—> 3<x u 

10.63 

12.30 

12.40 

14.21 

15.07 

12.70 

3 n u 

3<r g 

—¥ 2n u 

10.99 

12.30 

12.36 

13.04 

13.43 

12.71 

‘n u 

3<r g 

-¥ 2tt u 

10.98 

12.39 

12.44 

13.23 

13.45 

12.77 

X E+ 

3<7g 

—¥ 3cr u 

10.62 

12.43 

12.51 

14.31 

15.04 

12.82 

Ionization threshold: — chomo (eV) 




6.30 

14.94 


16.74 



MAD of excitation energies with respect to EOM-CCSD (eV) 

Valence 



0.48 

0.47 

0.35 

1.14 

1.65 

- 

Rydberg 



1.83 

0.34 

0.27 

1.17 

1.71 

- 

Total 



1.06 

0.41 

0.32 

1.15 

1.68 

- 


Maximum absolute deviation of excitation energies with respect to EOM-CCSD (eV) 


2.19 0.90 0.71 1.86 3.47 


Table II: Same as Table I for CO. 


State 

Transition 

TDKS 

TDRSH 

TDRSH+BSE2 

TDHF 

TDHF+BSE2 

EOM-CCSD 

Valence excitation energies (eV) 

3 n 

5ai(cr 

) —1 2ei(ir*) 

6.04 

6.10 

6.32 

5.85 

8.27 

6.45 

3 S + 

lei (tt 

) —» 2ei(-7r*) 

8.54 

8.45 

8.63 

7.79 

10.38 

8.42 


5ai(er 

) -> 2ei(-7r*) 

8.42 

8.68 

8.88 

9.08 

10.94 

8.76 

3 a 

lei (tt 

) 2ei(7r*) 

9.20 

9.13 

9.31 

8.74 

11.19 

9.39 

3 e~ 

lei (tt 

) -> 2ei( tt*) 

9.84 

9.80 

9.98 

9.73 

11.76 

9.97 


lei (tt 

) -> 2ei(7r*) 

9.84 

9.80 

9.98 

9.73 

11.82 

10.19 

X A 

lei (tt 

) —¥ 2ei(-7r*) 

10.33 

10.32 

10.50 

10.15 

12.05 

10.31 

3 n 

4ai (cr 

) -> 2ei(7r*) 

11.43 

11.96 

12.12 

13.31 

15.70 

12.49 

Rydberg excitation energies (eV) 


5ai(cr 

) ->• 6ai(cr) 

9.56 

10.34 

10.46 

11.18 

12.09 

10.60 

1 E+ 

5ai(cr 

) —1 6ai(cr) 

9.95 

11.12 

11.20 

12.27 

12.61 

11.15 

3 e+ 

5ai(cr 

) —1 7ai(<r) 

10.26 

11.08 

11.17 

12.42 

12.83 

11.42 

1 E+ 

5ai(cr 

) -> 7ai(<r) 

10.50 

11.30 

11.38 

12.79 

12.91 

11.64 

3 n 

5a^(a 

) 1 3ei(ir) 

10.39 

11.26 

11.34 

12.60 

13.20 

11.66 

x n 

5ai(cr 

) ->• 3ei(ir) 

10.50 

11.45 

11.52 

12.88 

13.21 

11.84 




Ionization threshold: — chomo 

levy- 






9.12 

13.49 


15.11 




MAD of excitation energies with respect to the EOM-CCSD calculation (eV) 


Valence 

0.33 

0.23 

0.16 

0.49 

2.02 

Rydberg 

1.19 

0.29 

0.22 

0.97 

1.42 

Total 

0.70 

0.26 

0.19 

0.69 

1.76 

Maximum absolute deviation of excitation 

energies 

with respect to EOM-CCSD (eV) 


1.34 

0.53 

0.36 

1.16 

3.22 


long-range response kernel, as well as the neglect of the 
GW2 quasiparticle correction for small enough value of 

M- 

B. Excitation energies of the N2, CO2, H2CO, and 
C2H4 molecules 

We now test the calculation of excitation energies with 
the complete proposed TDHF+BSE2 method, i.e. in¬ 
cluding now the short-range Hxc kernel in the linear- 
response part and neglecting the GW2 quasiparticle cor¬ 


rection. The excitation energies for each method and 
each molecule are given in Tables I-IV. Mean absolute de¬ 
viations (MAD) and maximum absolute deviations with 
respect to the EOM-CCSD reference are also given for 
valence, Rydberg, and all excitation energies. 

As already well known, TDKS with the LDA func¬ 
tional gives very underestimated Rydberg excitation en¬ 
ergies. TDRSH greatly improves the excitation energies 
for the Rydberg states and, to a lesser extent, for the va¬ 
lence states, resulting in total MADs of 0.41, 0.26, 0.13, 
and 0.14 eV for N 2 , CO 2 , H 2 CO, and C 2 H 4 , respectively. 
TDRSH also offers a more accurate description of valence 
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Table III: Same as Table I for H 2 CO. 


State 

Transition 

TDKS 

TDRSH 

TDRSH+BSE2 

TDHF 

TDHF+BSE2 EOM-CCSD 

Valence excitation energies (eV) 

j .4 2 

2^2 (n) 

-» 2bi (- 7 T* ) 

3.08 

3.17 

3.45 

3.76 

6.86 

3.56 

3 A 2 

2^2 (n) 

—V 26i (n *) 

3.70 

3.82 

4.11 

4.58 

7.37 

4.03 

3 Ai 

lbl ( 7 r) 

2 bi( 7 r*) 

6.35 

6.08 

6.39 

4.96 

8.30 

6.06 

3 b 1 

5 ai(cr) 

—y 26i (tt* ) 

7.77 

8.09 

8.40 

8.60 

12.28 

8.54 

Rydberg excitation energies (eV) 

3 B 2 

262(n) 

— > 6ai(cr) 

5.85 

6.83 

6.92 

8.17 

8.63 

6.83 

1 b 2 

262(11) 

—y 6ai(cr) 

5.93 

7.01 

7.08 

8.56 

8.72 

7.00 

3 b 2 

262(71) 

-* 7ai(er) 

6.96 

7.69 

7.81 

9.04 

9.85 

7.73 

3 A ± 

2 b 2 (n) 

—»■ 3 b 2 (<r) 

6.73 

7.77 

7.83 

9.24 

9.58 

7.87 

'b 2 

2b 2 (n) 

->■ 7oi(cr) 

7.04 

7.91 

8.00 

9.41 

9.78 

7.93 

1 A 1 

2b 2 (n) 

-> 3 62 (cr) 

6.78 

7.93 

7.97 

9.53 

10.01 

7.99 

1 a 2 

2b 2 (n) 

-> 36i(ir) 

7.55 

8.32 

8.39 

10.04 

10.26 

8.45 

3 a 2 

2b 2 (n) 

-> 36i(ir) 

7.58 

8.31 

8.38 

9.93 

11.07 

8.47 

3 b 2 

2b 2 (n) 

—y 8ai(cr) 

7.97 

8.90 

8.98 

10.21 

11.96 

8.97 

3 b 2 

2 b 2 (n) 

— y 8ai(cr) 

8.19 

9.17 

9.25 

10.86 

11.05 

9.27 




Ionization threshold: — chomo 

(eV) 






6.30 

10.33 


12.04 




MAD of excitation 

energies with respect to the EOM-CCSD calculation (eV) 


Valence 



0.47 

0.27 

0.17 

0.48 

3.15 

- 

Rydberg 



0.99 

0.07 

0.06 

1.45 

2.04 

- 

Total 



0.84 

0.13 

0.09 

1.17 

2.36 

- 

Maximum absolute deviation of excitation energies with respect to EOM-CCSD (eV) 




1.21 

0.45 

0.33 

1.59 

3.74 

- 




Table IV: Same 

as Table I for C 2 H 4 . 



State 

Transition 

TDKS 

TDRSH 

TDRSH+BSE2 

TDHF 

TDHF+BSE2 EOM-CCSD 

Valence excitation energies (eV) 

3 B lu 

1&3uM —» lb 2 g(7r* 

) 4.74 

4.35 

4.73 

3.54 

6.06 

4.41 

1 B 1 u 

lb 3u (7r) -*• lb 2 g(7r* 

) 7.91 

8.07 

8.38 

7.70 

9.11 

8.00 

3 B 1s 

lb 3 g(a) -y lb 2g (7T* 

) 7.18 

7.92 

8.04 

8.48 

10.43 

8.21 

1 Big 

lb 3g (cr) -> lb 2 g(7r* 

) 7.48 

8.04 

8.24 

9.23 

10.81 

8.58 

Rydberg excitation energies (eV) 

3 b 3u 

lf>3u(w) 4oi g (o-) 

6.59 

7.21 

7.35 

6.91 

7.37 

7.16 

1 B 3 u 

lb3u(w) ~y 4ai g (cr) 

6.65 

7.36 

7.48 

7.14 

7.43 

7.30 

3 Sig 

l& 3 u(^) —!> 2 b 2u (cr) 

6.98 

7.42 

7.78 

7.66 

8.10 

7.91 

3 52 g 

Ibsuf^) —¥ 3bi u (cr) 

7.10 

8.03 

8.11 

7.79 

8.07 

7.93 

1 Big 

lb 3u (ir) -> 2b 2u (cr) 

7.19 

7.92 

8.17 

7.75 

8.09 

7.97 

1 S 2g 

Ibsuf^) —¥ 3b iu (cr) 

7.15 

8.13 

8.20 

7.92 

8.07 

8.01 

3 Ag 

lb3u(w) -> 2b 3u (7r) 

8.03 

8.46 

8.60 

8.02 

8.64 

8.48 

X Ag 

lb 3 u(7r) -> 2b 3u (-7r) 

8.30 

8.87 

8.99 

8.61 

8.88 

8.78 

3 B 3u 

lf>3u(^) — ¥ 5oi g (cr) 

8.26 

8.97 

9.12 

8.74 

9.26 

9.00 

1 B 3 u 

lb 3 u (ir) — > 5oig(cr) 

8.28 

9.09 

9.20 

8.92 

9.13 

9.07 



Ionization threshold: — chomo 

JeV) 





6.89 

10.45 


10.23 



MAD of excitation energies with respect to the EOM-CCSD calculation (eV) 

Valence 


0.64 

0.24 

0.30 

0.52 

1.80 

- 

Rydberg 


0.71 

0.10 

0.17 

0.21 

0.14 

- 

Total 


0.69 

0.14 

0.21 

0.30 

0.62 

- 

Maximum absolute deviation of excitation energies with respect to EOM-CCSD (eV) 



1.10 

0.54 

0.38 

0.87 

2.23 

- 


and Rydberg excitation energies than TDHF. For a more 
intensive discussion of the performance of TDRSH, see 
Ref. 10. 

Both when starting from TDHF and TDRSH, the ad¬ 
dition of the BSE2 correlation correction always leads 
to larger excitation energies. In the case of TDHF, the 
BSE2 correction increases the valence excitation energies 
by about 2 or 3 eV, leading to largely overestimated va¬ 
lence excitation energies. For Rydberg states, the BSE2 
correction on top of TDHF is smaller (less than 1 eV) but 
also leads to systematically overestimated excitation en¬ 
ergies. Overall, TDHF+BSE2 considerably worsens the 
TDHF excitation energies. We thus conclude that, for 
these molecules, the relatively accurate results reported 


by Zhang et al. [52] crucially depend on using the GW2 
quasiparticle correction to the HF orbital energies. 

In the range-separated case, the long-range BSE2 cor¬ 
rection induces only a moderate increase of the va¬ 
lence excitation energies by about 0.2 to 0.4 eV, lead¬ 
ing for these states to MADs of 0.35, 0.16, 0.17, and 
0.30 eV for N 2 , C0 2 , H 2 CO, and C 2 H 4 , respectively. 
The TDRSH+BSE2 excitation energies of the Rydberg 
states are also systematically larger than the TDRSH 
ones by usually less than 0.1 eV, giving for these states 
MADs of 0.27, 0.22, 0.06, and 0.17 eV for N 2 , C0 2 , 
H 2 CO, and C 2 H 4 , respectively. Of course, the differ¬ 
ence in magnitude of the BSE2 correction in the range- 
separated case compared to the full-range case is to be 
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Figure 6 : Mean error versus standard deviation for the valence 
and Rydberg excitation energies of the N 2 , CO 2 , H 2 CO, and 
C 2 H 4 molecules calculated with linear-response TDKS (with 
the LDA functional), TDRSH and TDRSH+BSE2 (with the 
short-range LDA functional and n = 0.35 bohr -1 ), all within 
the TDA. The errors are calculated with respect to the EOM- 
CCSD excitation energies. The Sadlej+ basis set is used. 


mostly attributed to the substitution of the full-range 
two-electron integrals by the long-range ones. Since for 
the chosen value of the range-separation parameter /i 
of 0.35 bohr -1 , TDRSH mostly gives slightly underes¬ 
timated excitation energies of the considered systems, 
the long-range BSE2 correction overall slightly improves 
the excitation energies. More specifically, in comparison 
with TDRSH, TDRSH+BSE2 gives slightly smaller total 
MADs of 0.32, 0.19, and 0.09 eVfor N 2 , C0 2 , and H 2 CO, 
and a slightly larger MAD of 0.21 eV for C 2 H 4 . Also, for 
all the four molecules, TDRSH+BSE2 always gives the 
smallest maximum absolution deviation among all the 
methods, suggesting that TDRSH+BSE2 describes more 
reliably the excitation energies than the other methods. 

Finally, as a global summary of the results, Figure 6 
reports the mean error versus the standard deviation for 
the valence and Rydberg excitation energies of the four 
molecules for the different methods. For the valence exci¬ 
tation energies, going from TDKS to TDRSH mainly de¬ 
creases the standard deviation, while going from TDRSH 
to TDRSH+BSE2 decreases the mean error. For the Ry¬ 
dberg excitation energies, TDRSH provides a large im¬ 
provement over TDKS both in terms of mean error and 
standard deviation, while TDRSH+BSE2 gives a slightly 
smaller mean error than TDRSH. 


VI. CONCLUSION 

We have developed a range-separated linear-response 
TDDFT approach using a long-range frequency- 
dependent second-order Bethe-Salpeter correlation ker¬ 
nel. We have tested our approach using a perturba¬ 
tive resolution of the linear-response equations within the 
TDA for valence and Rydberg excitation energies of small 
atoms and molecules. The results show that the addition 


of the long-range correlation kernel overall slightly im¬ 
proves the excitation energies. 

More intensive tests should now be carried out with 
this long-range correlation kernel to better assess its per¬ 
formance. In particular, this long-range correlation ker¬ 
nel is expected to be appropriate for ( 1 ) calculating exci¬ 
tation energies of excited states with significant double¬ 
excitation contributions, ( 2 ) calculating charge-transfer 
excitation energies, and (3) calculating dispersion inter¬ 
actions in excited states. 

A number of further developments should also be ex¬ 
plored: adding the self-energy (or quasiparticle) contri¬ 
bution directly to the kernel, going beyond the TDA and 
the perturbative resolution of the linear-response equa¬ 
tions, and going beyond the second-order approximation. 
Finally, we note that the present work could be repeated 
using a linear decomposition of electron-electron inter¬ 
action [ 88 ], instead of a range separation, in order to 
construct a new TDDFT double-hybrid method which 
would be an alternative to the one commonly used based 
on CIS(D) [14], 
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Appendix A: Non-interacting four-point linear 
response function 

In this appendix, we derive the expression of the non¬ 
interacting linear-response function xo(a/,w) depending 
on two frequencies which is used in Eq. (32). 

The non-interacting four-point linear-response func¬ 
tion is defined in the time domain by 

Xo(xifi,x 2 f 2 ;xiti,x' 2 t 2 ) = 

-i G 0 (xiti, x 2 t' 2 )Go(x 2 f 2 , xi t[), (Al) 

where Go is the non-interacting one-electron Green func¬ 
tion. Using time-translation invariance and introduc¬ 
ing the time variables t\ = ti — t[, r 2 = t 2 — t 2 , and 
r = (ti +t [)/2 — (i 2 +t 2 )/2, Eq. (Al) becomes 

Xo(xi,x 2 ;x' 1 ,x 2 ;ti,t 2 ,t) = 

■n ( ' , Tl + T2 ( / T i + T i \ 

- *Go I xi,x 2 ,tH--—)G 0 (x 2 ,x 1 ,—-- t I , 

(A2) 

with G 0 (x 1 ,x' 1 ,ti — t[) = Go(x 1 t 1 ,x^). The triple 
Fourier transform of xo( x i, x 2 ; x^, x 2 ; ti, t 2 , t) is easily 


O 


A 

□ 


O Valence TDKS ■ A 

□ Valence TDRSH 
A Valence TDRSH+BSE2 
Rydberg TDKS 
■ Rydberg TDRSH 
A Rydberg TDRSH+BSE2 

.2 -1 -0.8 -0.6 -0.4 -0.2 0 

Mean error (eV) 
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found to be 


[1 /(<*;' - A) - l/(u/ - B)]/(A - B) gives 


Xo(xi,x 2 ;xi,x , 2 ;w , ,w' / ,a;) = 


dridr2dr e 


IUJ Ti ILJ T2 


e WT Xo (xi, x 2 ; xi, x 2 ; ti , r 2 , r) = 


-2iti5(d' - d") G 0 (xi,x' 2 ,u/ + ^ G 0 (x 2 ,xi,a/ - ^ 


(A3) 


where Go(xi, x' : , d) = J drie lUJTl Go(xi,x' 1 ,ri) is the 
Fourier transform of the Green function. The linear- 
response function xo(xi, x 2 ; x' 1; x 2 ; d', d) depending on 
two frequencies is then given by 

Xo(xi,x 2 ;x / 1 ,x 2 ;w , ,o;) 

= Xo(xi,x 2 ;xi,x 2 ;n = 0 ~,d',d) 

= Xo(xi,x 2 ;xi ,x' 2 ;uj',t 2 = 0 _ ,w) 

= J e * w "° + Xo( x i) x 2 ; x^ ,x' 2 ;d',d",d) 

= —i e lul ' 0+ G 0 (xi,x' 2 ,u/ + G 0 (x 2 ,xi,o J - ^ . 

(A4) 

Inserting in Eq. (A4) the Lehmann representation of the 
Green function, 


Xo( x i, x 2 ; x(,x 2 ; u/, d) = 


i e 


iu CP 


E 


i 


^(xi)^a(x i)y>a(x' 2 )</7j(x 2 ) 
id — (e a — £j) + z0+ 

1 


uj f -|- uj f 2 — £ a -|- i0 + uj r — uj /2 — £i — zO - * - 

Vl ( x i )<Pi (Xl )lfi* (x' 2 ) ip a (X 2 ) 


-V ■ 

-id - (e a - £i) + *0+ 
ia v 7 

l 

— d/2 — e a + *0+ d' + d/2 — £j — *0+ 


1 


E 


i 


y , a(x'i)yb(xi)^(^)y a (x 2 ) 
d - {e b - e a ) 

1 


UJ r + Cj/2 — £}) i0 + ci/ — cj/2 — H- 

Pi ( x i ) Pi (xi ) P* (x' 2 ) p* (x 2 ) 


E 


r /.. „/ V- Pa(xi)p:(xi) 

G 0 (x 1 ,x 1 ,a;)-^ a; _ £a+ . 0+ 


E 


Pi(xi)p*(x , i) 
d — £i — i 0 + ’ 


W - (£j - £i) 

1 


d' + d/2 — £j — i0 + d' — d/2 — £i — i0 + 


(A 6 ) 

(A 5 ) In Eq- (A 6 ), the first sum corresponds to the matrix ele¬ 
ment Xo ,ia,ia{<d',d) written in Eq. (32), while the second 
where i and a refer to occupied and virtual spin orbitals, sum corresponds to the matrix element xo,ai,ai(j-o' ,d) = 
respectively, and using the identity 1 /[(tj'—A)(u' — B)\ = Xo ,ia,ia(u', — d). 
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